Using ggplot

Please follow the instructions in the README file before starting this tutorial.

---

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

---

Install missing packages

• First check that we have all the packages that we need for the tutorial
• if any are missing then they will be installed

if (!require(tidyverse)) install.packages("tidyverse")
if (!require(readxl)) install.packages("readxl")
if (!require(RSQLite)) install.packages("RSQLite")
if (!require(DT)) install.packages("DT")
if (!require(plotly)) install.packages("plotly")
if (!require(ggsignif)) install.packages("ggsignif")
if (!require(ggdendro)) install.packages("ggdendro")
if (!require(patchwork)) install.packages("patchwork")

---

Loading the packages

• We will use the various packages of the tidyverse
• This includes the popular plotting library ggplot2
• Lets load the packages

library(tidyverse)

Loading the data

• We now need load some data to plot

• We can load this from most of the popular formats

• From a comma separated formated file (csv):

read_csv('data/stats.csv', show_col_types = FALSE)

• From a tab separated formated file (tsv):

read_tsv('data/stats.tsv', show_col_types = FALSE)

• From an excel file, in order to do this we need to first load the readxl package

library(readxl)

read_xlsx('data/stats.xlsx')

• We can also load data directly from a database, to access a local sqlite3 database, first load the RSQLite package (to load from MIcrosoft Access you could use the RODBC package)

library(RSQLite)

conn <- dbConnect(RSQLite::SQLite(), "data/stats.db")

dbGetQuery(conn, "SELECT * FROM stats")

• Above we loaded the data and printed it to the screen. But we need to keep the data in a variable:

stats <- read_xlsx('data/stats.xlsx')

• the data is now stored in the stats variable, we can verify this by just running the variable name

stats

• We can generate nicer looking tables for our reports using the DT package and the datatable function

library(DT)

datatable(stats)

• if we need to we can use the pipe symbol %>% to run other filtering and renaming functions on our data, here we use the filter() function to filter the data by the reads column. Notice that the number of rows reduces from 58 to 21

stats %>%
  filter(reads > 1000000)

• here we filter by both the reads and the coverage columns, we now only have 8 rows

stats %>%
  filter(reads > 1000000, coverage > 110)

• here we use the mutate() function to add a new column called total by calculating the sum of the coverage and the reads columns

stats %>%
  mutate(total = reads + coverage)

• here we use the rename() function to change the reads column name to counts

stats %>%
  rename(counts = reads)

---

Plotting with ggplot

Plotting a bar charts

• ggplot allows us to generate figures in layers, building them up piece by piece
• initially we generate the coordinate system
• ggplot() need the name of the variable containing the data and which variables to use in the plot

ggplot(stats, aes(x=isolate, y=reads)) 

• there is no actually plot here as we have not yet told ggplot which kind of plot we want
• we add the layers by simply using the ‘+’ operator. Here we specify a bar chart. By default geom_bar() will calculate the height of the bars from the data, but if we want to specify the data for the y-axis we need to supply stat="identity"

ggplot(stats, aes(x=isolate, y=reads)) +
  geom_bar(stat="identity")

• we can easily specifiy some colours

ggplot(stats, aes(x=isolate, y=reads)) +
  geom_bar(stat="identity", colour='black', fill='blue')

• common colours can be defined simply by name as above, but more control can be gained by specifying html web colours, a colour selector can be found here https://www.w3schools.com/colors/colors_picker.asp

ggplot(stats, aes(x=isolate, y=reads)) +
    geom_bar(stat="identity", colour='#ff0066', fill='#cc0000')

• We can change other stylistic parts of the chart with themes, we can set the general theme with some prebuilt styles e.g. theme_bw(), theme_minimal(), theme_classic(); more details and examples can be found at https://ggplot2.tidyverse.org/reference/ggtheme.html
• Lets try out the theme_bw() theme

ggplot(stats, aes(x=isolate, y=reads)) +
    geom_bar(stat="identity", colour='black', fill='blue') +
    theme_bw()

• the x-axis labels overlap and are unreadable, so lets rotate them using theme()

ggplot(stats, aes(x=isolate, y=reads)) +
  geom_bar(stat="identity", colour='black', fill='blue') +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle=-90)
  )

• the x-axis labels are not quite aligned to the tick marks, we can adjust them with the vjust and hjust parameters

ggplot(stats, aes(x=isolate, y=reads)) +
  geom_bar(stat="identity", colour='black', fill='blue') +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
  )

• we can now change the axes labels (labs) and add a title (ggtitle), change the font and font sizes of the axis labels

ggplot(stats, aes(x=isolate, y=reads)) +
  geom_bar(stat="identity", colour='black', fill='blue') +
  ggtitle("Read counts") +
  labs(x="Isolate name", y="Read count") +  theme_bw() +
  theme(
    axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5, size=10),
    axis.title  = element_text(size=15, family='Comic Sans MS')
  )

• if we prefer we can flip the whole chart

ggplot(stats, aes(x=isolate, y=reads)) +
  geom_bar(stat="identity", colour='black', fill='blue') +
  ggtitle("Read counts") +
  labs(x="Isolate name", y="Read count") +  theme_bw() +
  theme(
    axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
  ) +
  coord_flip()

• in the above examples we provided the stats data object to the ggplot() function. Very often its easier to pipe the data into the ggplot() function using %>%. As you can see, the stats variable name is now outside the ggplot() function

stats %>%
  ggplot(aes(x=isolate, y=reads)) +
    geom_bar(stat="identity", colour='black', fill='blue') +
    ggtitle("Read counts") +
    labs(x="Isolate name", y="Read count") +  theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    )

• This allows us to manipulate the data before we plot it, for example we can apply a filter and only plot certain samples. grepl() is just a pattern matching function, so filters isolates with the string SER in their name

stats %>%
  filter(grepl("SER",isolate)) %>%
  ggplot(aes(x=isolate, y=reads)) +
    geom_bar(stat="identity", colour='black', fill='blue') +
    ggtitle("Read counts") +
    labs(x="Isolate name", y="Read count") +
    theme_bw()

• now that we have a smaller plot, we can also add other useful information, like the actual values of the bars, using geom_text()

stats %>%
  filter(grepl("SER",isolate)) %>%
  ggplot(aes(x=isolate, y=reads)) +
    geom_bar(stat="identity", colour='black', fill='blue') +
    ggtitle("Read counts") +
    labs(x="Isolate name", y="Read count") +
    geom_text(aes(label=reads), vjust=-0.3, size=3.5) +
    theme_bw()

---

Plotting a scatter plot

• Lets plot a different type of graph, this time plotting reads against coverage

stats %>%
  ggplot(aes(x=reads, y=coverage)) +
    geom_point() +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    ) 

• we can easily add a smoothed regression line with geom_smooth()

stats %>%
  ggplot(aes(x=reads, y=coverage)) +
    geom_point() +
    geom_smooth() +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    ) 

`geom_smooth()` using method = 'loess' and formula 'y ~ x'

• if we need to we can limit the region of the plot to include, using xlim() and/or ylim()

stats %>%
  ggplot(aes(x=reads, y=coverage)) +
    geom_point() +
    geom_smooth() +
    xlim(c(0, 600000)) + 
    ylim(c(0, 100)) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    ) 

`geom_smooth()` using method = 'loess' and formula 'y ~ x'
Warning: Removed 21 rows containing non-finite values (stat_smooth).
Warning: Removed 21 rows containing missing values (geom_point).

• Finally if we wish we can make interactive versions of our plots using the plotly package
• try selection regions of the plot and zooming, using the menu ribbon on the top right of the plot

library(plotly)

coverage <- stats %>%
  ggplot(aes(x=reads, y=coverage)) +
    geom_point() +
    geom_smooth() +
    xlim(c(0, 600000)) + 
    ylim(c(0, 100)) +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    ) 

ggplotly(coverage)

`geom_smooth()` using method = 'loess' and formula 'y ~ x'
Warning: Removed 21 rows containing non-finite values (stat_smooth).

---

Exercise 1

• load the dataset stored in data/data.csv, plot a bar chart of county against popadults for the state of Ohio (OH)
• remember how to rotate the x axis labels to make them more readable,
• add a title to the plot and rename the axes titles

• Using the same dataset generate a scatter plot of area against poptotal
• Try adding a regression line using the glm method
• Try removing some of the outliers from the plot. Hint, use xlim() and/or ylim()
• Change the size and colour of the points? Hint: check the help with ?geom_line
• Colour the points by the state column
• Set the size of the points according to the popdensity column
• Move the legend position to the bottom of the chart. Hint: check the help with ?theme
• Add title and a subtitle
• Change the axes titles to ‘Population total’ and ‘County area’
• Improve the labels on the y-axis. Hint: look at ?scale_y_continuous

• Check out this site for lots of other example plots http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html

---

Plotting boxplots

• We will use one of the built in datasets for this section, the dataset shows the effect of Vitamin C on tooth growth via orange juice (OJ) or ascorbic acid (VC). Explore the data

ToothGrowth

• Lets do an initial scatter plot of dose against len colouring by supp

ToothGrowth %>%
  ggplot(aes(x=dose, y=len, color=supp)) +
    geom_point() +
    theme_bw()

• this is not that easy to visualise so lets jitter (separate) the points a little by replacing geom_point with geom_jitter

ToothGrowth %>%
  ggplot(aes(x=dose, y=len, group=supp, color=supp)) +
    geom_jitter(width=0.1) +
    theme_bw()

• In order to plot a boxplot the x-axis needs to be categorical data, so we can convert our dose values to characters using the mutate() and as.character() functions

ToothGrowth %>%
  mutate(dose=as.character(dose)) %>%
  ggplot(aes(x=dose, y=len, color=supp)) +
    geom_boxplot(color='black') +
    theme_bw()

• Sometimes it is useful to keep the individual points on the plot, we can just add back the jittered layer

ToothGrowth %>%
  mutate(dose=as.character(dose)) %>%
  ggplot(aes(x=dose, y=len, color=supp)) +
    geom_boxplot(color='black') +
    geom_jitter(width=0.1)  +
    theme_bw()

• we could even further break the boxplots down by the supp variable, by removing color="black" from geom_boxplot()

ToothGrowth %>%
  mutate(dose=as.character(dose)) %>%
  ggplot(aes(x=dose, y=len, color=supp)) +
    geom_boxplot() +
    geom_jitter(width=0.1)  +
    theme_bw()

• lets tidy up the plot, by labeling the axes and the legend properly

ToothGrowth %>%
  mutate(dose=as.character(dose)) %>%
  ggplot(aes(x=dose, y=len, color=supp)) +
    geom_boxplot(color='black') +
    geom_jitter(width=0.1)  +
    labs(x="Vitamin C dose", y="Tooth length") +
    theme_bw() +
    theme(
      legend.title = element_blank()    # this removes the legend title
    )

• And if we want to add some significance labels, we can make use of the ggsignif package

library(ggsignif)

ToothGrowth %>%
  mutate(dose=as.character(dose)) %>%
  ggplot(aes(x=dose, y=len, color=supp)) +
    geom_boxplot(color='black') +
    geom_jitter(width=0.1)  +
    labs(x="Vitamin C dose", y="Tooth length") +
    theme_bw() +
    theme(
      legend.title = element_blank()    # this removes the legend title
    ) +
    ylim(0,40) +
    geom_signif(comparisons = list(c("0.5", "1")), map_signif_level = TRUE, textsize = 6, y_position = 30, colour="black", annotation=c('**')) +
    geom_signif(comparisons = list(c("1", "2")), map_signif_level = TRUE, textsize = 6, y_position = 36, colour="black", annotation=c('*')) 

---

Plotting line plots

• Using the same tooth growth dataset, in order to draw a line plot we need to calculate the means and standard deviations of the data. Remember the data looks like this

ToothGrowth

• we can summarise the data by grouping the variables and calculating the mean/sd over these variables

summary <- ToothGrowth %>%
  group_by(supp, dose) %>%
  summarise(mean=mean(len), sd=sd(len))

`summarise()` has grouped output by 'supp'. You can override using the `.groups` argument.

summary

• it is now straightforward to generate a line plot

summary %>%
  ggplot(aes(x=dose, y=mean, color=supp)) +
    geom_line() +
    theme_bw()

• we can also use the calculated standard deviations to add error bars

summary %>%
  ggplot(aes(x=dose, y=mean, color=supp)) +
    geom_line() +
    geom_point() +
    geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.05) +
    theme_bw()

• the error bars overlap so we can actually shift (dodge) them slightly

summary %>%
  ggplot(aes(x=dose, y=mean, color=supp)) +
    geom_line() +
    geom_point() +
    geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.05,
                    position=position_dodge(0.05)) +
    theme_bw()

---

Exercise 2

• We will use another default dataset, mtcars. First take a look and explore

• Plot a boxplot with cyl on the x-axis and mpg on the y-axis
• include the raw data points, jittered and coloured by carb
• Tidy up the axes labels and remove the legend title

• Now plot a line graph showing the mean and the standard deviations of the mpg values, coloured by the number of carburetors

• Plot a similar plot comparing the number oc cylinders with the mpg but coloured by the number of gears

---

Plotting dendrograms

• We can easily perform hierarchical clustering on a dataset and then draw a tree
• here we use the built in dataset of USarrests, first explore the data

USArrests

• Now we can perform the clustering using the hclust function, this uses euclidean distances by default
• But the distance object (dist(USArrests)) can be any kind of similarity matrix, e.g pearson correlations

hc <- hclust(dist(USArrests), "ave")  # hierarchical clustering
hc

Call:
hclust(d = dist(USArrests), method = "ave")

Cluster method   : average 
Distance         : euclidean 
Number of objects: 50 

• Now we draw the tree using the ggdendrogram() function from the `ggdendro package

library(ggdendro)

p <- ggdendrogram(hc, rotate = TRUE, size = 2)
p

• remember that we can easily make the plot interactive with ggplotly()

ggplotly(p)

---

Combining plots into a figure

• In order to combine plots into a publication ready figure, lets plot both the previous reads and the coverage plots, only this time we will store the plots in the variables reads and coverage, notice that they won’t be printed when you run the code

reads <- stats %>%
  ggplot(aes(x=isolate, y=reads)) +
    geom_bar(stat="identity", colour='black', fill='blue') +
    ggtitle("Read counts") +
    labs(x="Isolate name", y="Read count") +  theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    )

coverage <- stats %>%
  ggplot(aes(x=isolate, y=coverage)) +
    geom_bar(stat="identity", colour='black', fill='#888888') +
    ggtitle("Coverage") +
    labs(x="Isolate name", y="Coverage") +
    theme_bw() +
    theme(
      axis.text.x = element_text(angle=-90, hjust=1, vjust=0.5)
    )

• but we can show the plots by just “running” the variable names

reads

coverage

• we can now make a combined plot usingt the patchwork package
• firstly side by side with the ‘+’ operator

library(patchwork)

reads + coverage

• if we want them above and below we change the ‘+’ operator to ‘/’

reads / coverage

• We can add subplot labels (A, B etc)

reads / coverage + plot_annotation(tag_levels = 'A') + plot_layout(guides = "collect")

• Note that we can include more plots in the figure. Here we generate four versions of the reads plot but changing the underlying theme. We can also add a new title, note that this is overriding the ggtitle() already in the reads plot. Can you spot the differences in the themes?

p1 <- reads + theme_bw() + ggtitle('theme_bw()')
p2 <- reads + theme_classic() + ggtitle('theme_classic()')
p3 <- reads + theme_minimal() + ggtitle('theme_minimal()')
p4 <- reads + theme_dark() + ggtitle('theme_dark()')

( p1 + p2 ) / ( p3 + p4 ) + plot_annotation(tag_levels = 'A')

---

Exercise 3

• Generate a combined figure including a boxplot and a lineplot comparing the number of cylinders with the mpg as generated above, label the subplots using a, b

---

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

---

Session details

• Generate document version details

sessionInfo()

R version 4.1.1 (2021-08-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS:   /opt/R/4.1.1/lib64/R/lib/libRblas.so
LAPACK: /opt/R/4.1.1/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8       
 [4] LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggsignif_0.6.3  ggdendro_0.1.23 gapminder_0.3.0 gganimate_1.0.7 plotly_4.10.0   patchwork_1.1.1
 [7] DT_0.23         readxl_1.3.1    forcats_0.5.1   stringr_1.4.0   dplyr_1.0.7     purrr_0.3.4    
[13] readr_2.0.1     tidyr_1.1.3     tibble_3.1.4    ggplot2_3.3.5   tidyverse_1.3.1 RSQLite_2.2.14 

loaded via a namespace (and not attached):
 [1] nlme_3.1-152      fs_1.5.0          lubridate_1.7.10  bit64_4.0.5       progress_1.2.2   
 [6] httr_1.4.2        tools_4.1.1       backports_1.2.1   bslib_0.3.1       utf8_1.2.2       
[11] R6_2.5.1          DBI_1.1.1         lazyeval_0.2.2    mgcv_1.8-36       colorspace_2.0-2 
[16] withr_2.4.2       tidyselect_1.1.1  prettyunits_1.1.1 bit_4.0.4         compiler_4.1.1   
[21] cli_3.0.1         rvest_1.0.1       xml2_1.3.2        labeling_0.4.2    sass_0.4.1       
[26] scales_1.1.1      digest_0.6.27     rmarkdown_2.11    pkgconfig_2.0.3   htmltools_0.5.2  
[31] dbplyr_2.1.1      fastmap_1.1.0     htmlwidgets_1.5.4 rlang_0.4.11      rstudioapi_0.13  
[36] jquerylib_0.1.4   farver_2.1.0      generics_0.1.0    jsonlite_1.7.2    crosstalk_1.2.0  
[41] vroom_1.5.5       magrittr_2.0.1    Matrix_1.3-4      Rcpp_1.0.7        munsell_0.5.0    
[46] fansi_0.5.0       lifecycle_1.0.0   stringi_1.7.4     yaml_2.2.1        MASS_7.3-54      
[51] plyr_1.8.7        grid_4.1.1        blob_1.2.2        parallel_4.1.1    crayon_1.4.1     
[56] lattice_0.20-44   haven_2.4.3       splines_4.1.1     hms_1.1.0         knitr_1.34       
[61] pillar_1.6.2      reprex_2.0.1      glue_1.4.2        evaluate_0.14     data.table_1.14.0
[66] modelr_0.1.8      vctrs_0.3.8       tzdb_0.1.2        tweenr_1.0.2      cellranger_1.1.0 
[71] gtable_0.3.0      assertthat_0.2.1  cachem_1.0.6      xfun_0.26         broom_0.7.9      
[76] viridisLite_0.4.0 memoise_2.0.1     ellipsis_0.3.2